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ABSTRACT 


— *  The  recording  and  subsequent  analysis  of  evoked  potential 
activity  has  proven  useful  for  the  evaluation  of  neural  dysfunction 
resulting  from  impact  acceleration  injury  involving  the  head 
and  neck.  In  animal  impact  acceleration  experiments  involving 
Rhesus  monkeys,  somatosensory  evoked  potentials  showed  an  Increase 
in  latency  following  nonlethal  experiments.  In  order  to  assess 
quantitatively  and  objectively  the  amplitude  and  duration  of 
the  latency  effect  following  impact  acceleration,  a  nonlinear 
mathematical  model  has  been  proposed.  This  technical  report 
describes  a  nonlinear  regression  procedure  for  fitting  the 
proposed  model  directly  to  empirical  latency  data.  A  FORTRAN 
computer  program  listing  is  provided.  CT — ' — 
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1.  INTRODUCTION 


An  important  consequence  of  impact  injury  involving  the  head 
and  neck  is  the  disruption  of  normal  functioning  of  the  central 
nervous  system.  As  part  of  its  research  effort  on  impact  accel¬ 
eration  injury  prevention,  the  Naval  Biodynamics  Laboratory  (NBDL) 
has  conducted  a  series  of  experiments  designed  to  test  the 
neurophysiological  effects  of  indirect  or  inertial  head-neck  impact 
acceleration.  In  these  experiments,  unanesthetized  Rhesus  monkeys 
were  subjected  to  peak  sled  accelerations  in  the  -X  direction  and 
somatosensory  evoked  potentials  (EPs)  were  recorded  prior,  during, 
and  subsequent  to  impact.  These  experiments  are  discussed  in 
detail  in  Berger  and  Weiss  [  1] . 

A  primary  objective  of  these  experiments  was  to  determine 
the  extent  to  which  impact  produced  shifts  in  latency  of  various 
peaks  of  the  (cervical)  evoked  potentials.  Shifts  in  latency  of 
each  peak  were  plotted  as  a  function  of  time  over  the  experiment, 
relative  to  impact.  Figure  1  gives  an  example  of  such  a  plot. 

In  order  to  quantitatively  assess  the  amplitude  and  duration 
of  the  latency  effect  following  impact,  the  following  exponential 
model  was  proposed: 

y  *  B  +  St  +  h(t)D  +  h(t)Aexp(t/T)  +  £(t)  (1) 

where : 

y  is  the  value  of  the  shift  in  latency  with  respect  to 
the  preimpact  baseline  average  evoked  potential  (AEP) 


CERVICAL  AEP  LATENCY  SHIFT 

RUN  LXSOOt.  7*6  M/81.  P1AK  6.46  MS.  N-10 


Figure  1.  Example  Plot  of  the  Shift  in  Latency 
as  a  Function  of  Time,  Relative  to  Impact.  (Taken 
from  Berger  and  Weiss  [1].) 


>  s 


t  is  the  time  relative  to  impact 

S  is  the  slope  of  the  linear  baseline  shift 

B  is  the  amplitude  value  of  the  baseline  instantaneously 
prior  to  impact 

D  is  the  "permanent"  shift  induced  by  impact 

A  is  the  exponential  amplitude 

T  is  the  exponential  time  coefficient 

h(t)  is  the  unit  step  function  with  value  zero  for  t<0, 
and  unity  for  t>0 

e(t)  is  the  noise. 


In  Berger  and  Weiss  [1],  the  linear  baseline  shift 
parameters,  S  and  B,  were  determined  by  applying  simple  linear 
regression  to  the  preimpact  data.  The  linear  equation  was 
extrapolated  in  the  postimpact  region,  and  residual  values  were 
computed  by  subtracting  out  the  linear  equation  component.  However, 
a  regression  procedure  which  would  fit  the  model  D  +  Aexp(t/T) 
to  the  postimpact  residuals  was  not  available.  Consequently, 
an  additional  exponential  term  was  introduced  in  lieu  of  the  D 
term,  and  the  residuals  were  then  subjected  to  a  polyexponential 
regression  procedure.  More  specifically,  the  (postimpact)  residuals 
were  fit  to  a  regression  function  of  the  following  form: 


A^exp(t/Tj )  +  A2exp(t/T2). 

The  estimate  of  the  "permanent"  shift  due  to  impact,  D,  was 

-3- 


obtained  from  the  amplitude  coefficient  corresponding  to  the 
exponential  term  having  a  time  coefficient  (T^)  that  was  very 
long  in  relation  to  the  duration  of  the  analysis  interval.  In 
some  cases  a  constant  value  was  added  to  the  data  to  ensure  that 
term  with  a  long  time  coefficient  would  be  obtained. 

In  this  Desmatics  technical  report  we  present  an 
alternative  computational  procedure  for  fitting  model  (1) 
directly  to  empirical  latency  data.  In  addition,  we  give  a 
FORTRAN  program  implementing  the  proposed  computational  method. 
Finally,  we  include  an  example  which  illustrates  a  direct 
application  of  the  methodology. 


2.  METHODS 


We  can  rewrite  model  (1)  as 

y±  ■  f^,  0)  +  e(t±)  (i  -  1,  2 . n) 

where  y^  is  the  measured  latency  shift  corresponding  to  the 
known  time  t^,  e(t^)  is  the  random  error  component,  n  denotes 
the  total  number  of  data  points,  and  the  vector  0  -  (S,  B,  D, 

A,  T)  consists  of  the  5  parameters  to  be  estimated.  Assuming 
that  the  errors  e(t^)  all  have  zero  mean,  the  expected  responses, 
denoted  by  Efy^),  can  be  ’-epresented  as: 

E(y±)  “  f(ti’ 

where  f  is  of  the  form 

f(t,  0)  -  B  +  St  +  h(t)D  +  h(t)Aexp(t/T) .  (2) 

If  we  assume  further  that  the  errors  e(t^)  all  have  the  same 
variance  and  are  uncorrelated,  then  it  is  reasonable  to  use  the 
least  squares  method  to  fit  equation  (2)  to  the  data.  That  is, 

A 

we  wish  to  find  the  estimate  of  0,  say  0,  which  minimizes  the  sum 
of  squares 

n  2 

Q(0)  -  T.  [y  -  f(t.,  0)T. 
i-1 
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The  proposed  procedure  is  based  on  the  following 
approach.  For  known  values  of  S  and  B,  say  s  and  b,  we  define 
y^Cs.b)  “  y^  -  b  -  st^.  Note  that 


E[yi(s,b)] 


0 

D  +  Aexp(t^/T) 


if  t±  <  0 
if  t^  >  0. 


We  will  assume  henceforth  that  the  data  has  been  ordered  and 

indexed  so  that  t,  <  t»  <  t_<  ...<  t  .  In  addition,  we  let  n, 

1  4  J  n  1 

and  n^  denote  the  total  number  of  preimpact  (t^  <  0)  and  postimpact 
(t^  >  0)  data  points,  respectively. 

It  is  clear  that  for  known  values  of  S  and  B  the 
parameters  D,  A,  and  T  can  be  estimated  from  an  analysis  of  the 

/v 

postimpact  residuals  y^Cs.b),  i  >  n^.  Accordingly,  let  D(s,b), 

A  A 

A(s,b),  and  T(s,b)  denote  the  least  squares  estimates  obtained 
from  fitting  the  exponential  function  D  +  Aexp(t^/T)  to  the 
residuals  y^(s,b),  i  >  n^.  [We  assume  temporarily  that  we  have 

A  A  A 

a  method  at  our  disposal  for  computing  D(s,b),  A(s,b),  and  T(s,b).] 
If  we  define 


V(s,b)  -  Q(8(s,b))  «  Z  [y,  -  f  (t, ,6(s,b))  ]2 

i  =  l  1 

A  AAA 

where  0(s,b)  ■  (s,  b,  D(s,b),  A(s,b),  T(s,b))  we  have  the  following 
key  result: 
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Let  s  and  b  denote  the  values  of  S  and  B,  ^ 
respectively,  which  minimize  V(s,b).  Then  0, 
the^value  of  0  which  minimizes  Q(0),  is  given 

by  0  -  (s,  b,  D(s,b) ,  A(s,b),  T(s,b>). 

The  significance  of  this  result  is  that  the  stated  problem  of 
minimizing  a  five-variable  function  can  be  reduced  to  that  of 
minimizing  the  two-variable  function  V(s,b).  Of  course,  an 
analysis  of  V  will  require  that  we  have  available  a  method  for 
fitting  the  exponential  model  D  +  Aexp(t^/T)  to  the  residuals 
y^(s,b)  for  some  s  and  b,  where  i  >  n^.  In  Section  3  a  general 
method  for  fitting  an  exponential  growth  and  decay  model  with  an 
asymptote  will  be  described.  In  Section  4  a  method  for 
minimizing  V(s,b)  will  be  presented. 


3.  EXPONENTIAL  REGRESSION  WITH  AN  ASYMPTOTE 


In  this  section  we  address  the  problem  of  finding  the 

A  AAA 

estimate  of  X,  say  X  =  (X^,  X^.  X^)»  which  minimizes  the  sum  of 
squares 

F(X)  -  E  [z  -  w(t  ,X)]2,  (3) 

i>nt 

where  wCt^.X)  =  X^  +  X^expCX^t^)  and  =  y^Cs.b),  with  t^  >  0. 

The  summation  in  (3)  contains  n^  terms,  one  for  each  postimpact 
data  point.  We  will  use  a  modified  Newton's  method  to  compute 

A 

the  vector  X. 

Newton's  method  is  a  widely  used,  iterative  method  for 
minimization  and  requires  use  of  both  the  gradient  vector  and  the 
Hessian  matrix  in  computations.  The  gradient  vector  consists  of 
all  first-order  partial  derivatives  of  the  objective  function  with 
respect  to  the  unknown  parameters;  the  Hessian  matrix  consists 
of  all  second-order  partial  derivatives.  [For  the  sake  of  discussion, 
we  will  assume  that  all  first-  and  second-order  derivatives  exist 
and  are  continuous.] 

To  provide  further  description  of  Newton's  method,  consider 
a  Taylor  series  expansion  of  the  objective  function,  F(X) ,  about 
the  point  X^.  This  takes  the  form 

F(X)  =  F^)  +  (X  -  Xi)'g1  +  Jj(X  -  X^'G^X  -  \±)  +  Rt 
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approximate  the  function  V(s,b)  in  the  design  region.  If  the 
fitted  surface  is  an  adequate  approximation  of  V,  then  analysis  of 
the  fitted  surface  will  be  approximately  equivalent  to  an  analysis 
of  the  actual  function.  An  indication  of  lack  of  fit  about  the 
fitted  surface  can  be  obtained  by  considering  variation  between 
the  observed  responses  (i.e.,  function  evaluations)  and  the  predicted 
responses  (i.e.,  fitted  values)  at  the  nine  points  which  comprise 

A  A  A  A 

the  CCD.  The  predicted  responses  v  =  (v^,  v^.-.-Vg)'  are  given 
by 


/s 

V  -  X§, 

A 

where  X  and  6  are  as  defined  previously.  A  statistic  which  can  be 
used  to  measure  lack  of  fit  is  the  mean  absolute  residual  (MAR) 
defined  by 

9 

MAR  »  (  L  I v.  -  v  | )  /9 
i=l 

where  v^ ,  v^.-.-.v^  are  the  nine  observed  responses.  The  deviations 
v^  -  v^  represent  the  difference  between  the  observed  and  fitted 
values  based  on  the  second-order  regression.  The  closer  are  the  v^ 

A 

to  the  v^,  the  smaller  is  the  value  of  MAR,  and  hence  the  better  the 
approximation  by  the  second-order  response  surface. 
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minimum  on  the  estimated  response  surface  is  given  by 


A 


X 


C 


A  A 


where  C 


(8j.  82)'  and 


B5/2 


B5/2 


In  the  original  space,  we  have 


s  =*  s  +  F  x, 
s  1 


and 

/\  •*« 

b  -  b  +  F,x„. 

b  Z 

It  should  be  noted  that  under  proper  scaling  of  the  design 

A  A 

region,  the  point  (x^,  x^)  should  be  located  somewhere  within  the 

~2  ~2 

design  region  (i.e.,  x^  +  <  2).  However,  if  this  point  should 

fall  outside  the  design  region,  then  a  new  CCD  design  should  be 

A  A 

set  up  with  (s,  b) ,  as  defined  above,  as  the  center. 

Lack  of  Fit 


The  second-order  response  surface  given  in  (A)  is  used  to 
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is  fit  to  the  nine  responses  (i.e.,  function  evaluations)  v^» 

Vg  corresponding  to  the  nine  points  in  the  CCD.  In  matrix  notation 
we  have 


v  -  X6, 


where  v  -  (v^  v?> . . .  ,v9> '  ,§  -  (3Q. 


, 85) ' .  and 


X  - 


1  1 
1  -1 
1  -1 
1  1 
1  0 
1  0 
1  ~/T 

1  /r 
1  0 


1  1 

1  1 

-1  1 

-1  1 

0  0 

JT  0 

0  2 

0  2 

-/T  0 


1  1 

1  -1 

1  1 

1  -1 

0  0 

2  0 

0  0 

0  0 

2  0 


It  follows  directly  that  the  estimate  of  the  regression  coefficients 
is  given  by 


b  -  (x’xrVv. 


We  can  also  show  that  the  point  x  ■  (x^,  x^)  which  results  in  the 
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where  F  and  F_  represent  scaling  factors.  In  terms  of  the  coded 
s  b 

variables,  we  have 

8  «  I  +  F  X. 

s  1 

and 


b  -  b  +  F,  x  . 

b  2 

Thus,  for  example,  the  coded  design  point  (x^,  X£>  ■  (/27  0)  corres¬ 
ponds  to  the  uncoded  design  point  (s,b)  ■  (I  +  / 2F  ,b) . 

s 

Initial  parameter  estimates  of  the  linear  baseline  shift 
parameters,  S  and  B,  can  readily  be  determined  by  applying  simple 
linear  regression  to  the  preimpact  data.  From  this  linear  regression 
one  can  also  determine  the  standard  errors  of  the  estimates  i  and  b. 
These  standard  errors  can  be  used  to  scale  the  design  region.  It 
seems  reasonable  to  define  the  scaling  factors  F  and  F.  as  one- 

S  D 

tenth  the  respective  standard  errors  of  s  and  b.  Additional  comments 
on  the  scaling  of  the  design  region  will  be  given  later  in  this 
section. 

As  we  indicated  earlier,  a  second-order  response  surface  in  the 
coded  variables 


vi  "  B0  +  6lXli  +  62X2i  +  ^3Xli  +  64X2i  +  B5XliX2i 
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explores  the  parameter  space  by  means  of  a  geometric  configuration 
of  points.  Function  evaluations  are  made  at  each  point  in  the 
design  configuration,  and  a  second-order  response  surface  is  fit 
to  the  responses  (i.e.,  evaluations)  to  approximate  V(s,b)  in  the 
design  region.  An  analysis  of  the  fitted  surface  is  then  made 
to  determine  the  point  which  results  in  the  minimum  of  the 
estimated  response  function. 

/v  A 

The  statistical  design  we  will  use  to  determine  (s,b)  is 
called  a  central  composite  design  (CCD).  A  CCD  for  two  variables, 
x^  and  x^,  is  illustrated  geometrically  in  Figure  3.  The  design 
consists  of  the  nine  following  coded  design  points:  (1,1),  (-1,1), 
(-1,-1),  (1,-1),  (0,0),  (0,ei),  (-ot,0),  (cx,0) ,  and  (0,-a),  where 
a  *•  JT7  The  first  four  points  are  the  usual  factorial  points  for 
fitting  a  first  order  model  for  two  variables.  The  fifth  point  is 
the  center  point  of  the  design,  and  the  four  remaining  points  are 
the  axial  points  of  the  square.  Notice  that  each  variable  is 
measured  at  five  coded  levels:  0,  1,  -1,  a,  and  -a. 

The  center  point  of  the  CCD  design  corresponds  to  our  initial 
estimates,  which  we  will  denote  by  2  and  b,  of  S  and  B.  The  coded 
variables  x^  and  x2  are  defined  as 

x  -  (s  -  s)/F 
1  s 

and 

x2  -  (b  -  b)/Ffe, 
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In  this  section  we  consider  the  problem  of  estimating  the  para 
meters  S  and  B,  the  slope  of  the  linear  baseline  shift  and  the 
amplitude  value  of  the  baseline  instantaneously  prior  to  impact, 
respectively.  In  Section  3  we  assumed  fixed  known  values,  s  and  b, 
for  these  parameters  and  outlined  an  iterative  method  for  obtaining 
least  squares  estimates  of  the  postimpact  parameters  D,  A,  and  T, 
based  on  an  analysis  of  the  residuals  y^(s,b)  ■  y  -  b  -  st^,  with 

A  A 

t  >  0.  Following  the  notation  used  in  Section  2,  let  D(s,b),  A(s,b) 

A 

and  T(s,b)  denote  the  least  squares  estimates  obtained  from  fitting 
the  exponential  function  D  +  AexpCt^/T)  to  the  residuals  y^Cs.b), 

A  A 

i>n^.  We  wish  now  to  find  the  estimates  S  and  b  which  minimize 
the  function 

V(s,b)  -  Q(0(s,b))  -  [yi  -  f^.QCs.b))]2 

A  A  A  A 

where  6(s,b)  «  (s,  b,  D(s,b),  A(s,b),  T(s,b))  and  the  summation 
extends  over  all  the  data. 

Because  V(s,b)  cannot  be  differentiated  analytically  and 
because  function  evaluations  (of  V)  are  computationally  demanding, 
to  determine  s  and  £  we  will  use  a  numerical  method  which  (1)  makes 
no  use  of  any  of  the  derivatives  of  V,  and  (2)  is  highly  efficient 
with  regard  to  the  number  of  function  evaluations  it  requires.  The 
proposed  method  is  based  on  the  use  of  a  statistical  design  which 
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interval  in  which  the  minimum  is  known  to  lie  by  a  factor  of  r. 

The  width  of  the  target  interval  at  each  iteration  is 
called  the  tolerance.  A  prescribed  accuracy,  a,  can  be  achieved 
for  estimating  A^  by  repeating  the  golden  section  procedure  until 
the  tolerance  is  less  than  or  equal  to  a. 


interval  [x0,  x  ].  A  function  f(x)  is  unimodal  in  the  interval 
X/  u 

* 

[x0 ,  x  ]  if  it  has  a  unique  minimum  (maximum) ,  occuring  at  x  , 

Aj  U 

in  the  interval,  is  a  decreasing  (increasing)  function  in  the 

it 

interval  [x^,  x  ],  and  is  an  increasing  (decreasing)  function  in 

A 

the  interval  [x  ,  x  ]. 

u 

A  complete  discussion  of  the  method  of  golden  section  can 
be  found  in  Keifer  [3].  However,  to  illustrate  the  basic  technique, 
consider  a  unimodal  function  f(u)  with  a  minimum  in  the  interval 
[Uj,  u^]  and  consider  the  four  points  u^  <  <  u^  <  u^  that 

satisfy  the  spacing  equations 


“S’  U1 


r(u4  -  uL) 


where  r  -  2/(1  +  =  .618034.  The  two  limits  u^  and  u^  are 

fixed  in  advance  and  are  known  to  bracket  the  minimum. 

By  testing  the  relative  values  of  f(u^),  fO^),  f(u^),  and 
f(u4),  it  is  possible,  since  f  is  unimodal  [see  Figure  2,  for 
example],  to  determine  in  which  of  the  two  intervals  [u^,  u^] 
or  [Ujt  u4]  the  minimum  lies.  Without  loss  of  generality  (note 
that  the  Intervals  are  of  Identical  length) ,  assume  that  the 
minimum  lies  in  the  interval  (u^,  u^).  The  process  is  then 
repeated  in  the  target  interval  (u^,  u^)  by  considering  the 
addition  of  a  new  interior  point  which  satisfies  the  basic 
spacing  equations.  Hence  each  subsequent  iteration,  which 
involves  one  function  evaluation,  reduces  the  previous  target 
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convergence  process,  especially  with  quadratically  convergent 
methods  such  as  Greenstadt's.  We  will  now  proceed  to  describe 
a  general  method  for  obtaining  initial  parameter  estimates  of 
X^ ,  X ^ ,  and 

Once  again,  consider  the  regression  function 
w^.X)  -  Xj  +  X2exp(X3t1), 

which  we  wish  to  fit  to  the  set  of  data  z ^  -  y^s.b),  i  >  n^. 

% 

Observe  that  for  a  fixed  known  value  of  Xy  say  Xy  estimates 

r\j  r\j 

of  Xj  and  Xy  call  them  X^  and  Xy  can  easily  be  determined  via 

'Xj  % 

a  simple  linear  regression  of  the  z ^  on  Xj  +  X^^  where  t^  ■ 
expCX^t^).  Accordingly,  define 

<\j  m  f\,  %  *\,  2 

W(X3)  -  l  [z±  -  Xj  -  X2exp(X3ti>]  , 
i>n^ 

where  X^  and  X2  depend  on  Xy  Thus,  an  initial  estimate  of  X3 
can  be  obtained  by  finding  the  value  of  X3  which  minimizes  W$3). 

%  'Vi 

The  corresponding  values  of  X^  and  X2  will  provide  initial 
parameter  estimates  of  X^  and  Xy 

To  minimize  the  function  W(X3),  we  will  use  the  method  of 
golden  section,  which  is  a  commonly  used  method  for  solving 
optimization  problems  in  one  variable.  This  method  assumes  that 
the  objective  function  is  a  unimodal  function  in  a  specified 


d(i,j) 


!1  for  j  -  1 

exp(X^t^)  for  j  -  2 

ti^2exP(^3ti>  for  J  "  3 


and 

!tiexP(X3ti) 
t^X2exp(X3ti) 

° 


for  j  «  2,  k  *  3,  and  j  *  3, 
for  j  ■  3,  k  ■  3 
otherwise. 


Using  these  results,  the  components  of  the  gradient  vector  of  F(X) 
have  the  form 


3F(X)/9X  -  -2  l  (z  -  w(t.,X))d(i,j), 

2  i>nL  1 

and  the  components  of  the  Hessian  matrix  of  F(X)  have  the  form 


32F(X)/3X  3X  «  -2  Z  (z  -  w(t  ,X))h(i,j,k)  +  2  Z  d(i,  j)d(i,k) . 
J  i>nL  1  1  i>nx 

Initial  Parameter  Estimates 


As  we  have  already  noted,  the  proposed  method  of  estimating 
X  requires  initial  parameter  estimates  of  X^,  X2»  and  X^. 

Reasonably  good  initial  estimates  should  be  provided,  however,  since 
poor  starting  estimates  can  have  a  detrimental  effect  on  the 


eigenvector,  with  v^v^  “  1*  Define  the  matrix  as 


where  -  max  (  |y^|,  <$)  and  6  is  a  small  positive  number. 

h 

The  matrix  is  positive  definite  with  inverse  given  by 


W*)'1  -  ^ 


Under  the  proposed  modification,  the  general  iterative  scheme 
is  given  by 


-i+1  “  -i  "  °‘i(-i)  §i 

where  the  scalar  a^,  a^>0,  is  a  step  size  factor  introduced  to 
help  avoid  stepping  across  a  minimum.  In  general,  is  selected 
so  that  ^(^:£+j)<  F(X^);  usually,  =  1  will  suffice.  Iterations 
are  carried  out  according  to  the  general  scheme  until  none  of 

A 

the  components  of  the  vector  X  change  by  more  than  a  negligible 
amount. 

For  use  in  G  eenstadt's  method,  we  will  now  derive  the 

gradient  vector  and  Hessian  matrix  of  the  objective  function  F(X) 

given  in  equation  (3).  Let  d(i,j)  -  3w(t^,X)/3Xj  and  h(i,j,k)  * 

2 

3  w(t^,X) /3X^3Xj .  It  can  easily  be  checked  that 
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I 


p 


where  g^  is  the  gradient  vector  evaluated  at  X^,  is  the  Hessian 
matrix  evaluated  at  X^,  and  R^  is  the  remainder.  The  terms  in 
the  Taylor  expansion,  ignoring  R^,  define  a  quadratic  function 
which  will  approximate  F(X)  when  X  is  close  to  X^.  The  stationary 
point  (i.e.,  a  point  at  which  the  gradient  vector  vanishes)  of 
this  quadratic  function  is  given  by  the  solution  to  the  linear 
system  of  equations 


2^^  ”  "  “Si* 

If  is  positive  definite  and  therefore  invertible,  this  suggests 
the  general  iterative  scheme 

-i+1  "  -i  '  -i  Si* 

This  recursive  relation  defines  Newton's  method,  although  one 
must  generally  provide  an  initial  starting  point  X^.  The  method 
normally  displays  a  quadratic  rate  of  convergence  near  a  stationary 
point. 

A  potential  problem  with  Newton's  method  is  that  the  Hessian 


matrix  may  not  be  positive  definite  at  each  iteration;  consequently, 
several  modifications  have  been  suggested.  Greenstadt  [2]  suggested 


the  following  variant  of  Newton's  method.  Let  y^  be  the  jth 


eigenvalue  of  the  p  x  p  matrix  G^  and  v  its  corresponding 
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— * "*  „*  **  • 0 


Jr  ^ 


5.  AN  EXAMPLE 


In  this  section  we  illustrate  by  means  of  an  example  a 
direct  application  of  the  methodology  developed  in  Sections  2,  3, 
and  4.  In  the  Appendix  we  give  a  FORTRAN-77  listing  of  the 
proposed  nonlinear  regression  method.  This  program  runs  interactively 
and  was  implemented  on  the  Zenith-100  microcomputer.  The  program 
can  easily  be  converted  to  run  on  a  mainframe  computer  or  another 
microcomputer  which  supports  a  FORTRAN  compiler. 

Cervical  EP  latency  data  was  obtained  from  the  Texas  Research 

Institute  of  Mental  Sciences.*  Latency  data  of  a  selected  EP  peak 

2 

from  Experiment  #LX3009  (796  m/sec  )  is  presented  in  Table  1.  For 
this  data,  each  test  AEP  spanned  10.2627  seconds  and  was  computed 
by  averaging  a  sequence  of  50  EPs  with  a  stimulus  sampling  rate  of 
4.872  Hz.  Latency  values  are  given  starting  at  six  minutes  prior 
to  impact,  continuing  through  impact,  and  ending  fifteen  minutes 
after  impact.  The  data  is  plotted  in  Figure  4.  It  should  be  noted 
that  the  three  data  points  marked  by  an  asterisk  (*)  in  Table  1 
were  Identified  as  spurious  data  values  and  omitted  from  subsequent 
analyses. 

A  baseline  latency  value  of  3.865  ms  was  computed  by  averaging 
the  preimpact  latency  values.  The  shift  in  latency  at  each  time 
point  was  determined  by  calculating  the  difference  between  the 

*The  author  is  indebted  to  Mr.  William  D.  Burton,  Systems  Analyst  at 
the  Texas  Research  Institue  of  Mental  Sciences,  for  supplying  the 
latency  data. 
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Relative 

Time 

(sec) 

Latency 

(ms) 

Latency 

Shift 

(ys) 

Relative 

Time 

(sec) 

Latency 

(ms) 

Latency 

Shift 

(ys) 

Relative 

Time 

(sec) 

Latency 

(ms) 

Latency 

Shift 

(ys) 

-356. 000 

3.  875 

10. 

64. 772 

4.025 

160. 

485. 543 

3.925 

60. 

-345. 737 

3.900 

35. 

75. 034 

3.950 

85. 

495. 806 

3.900 

35. 

-335. 475 

3.875 

10. 

85. 297 

3.  950 

85. 

506. 069 

3.900 

35. 

-325.212 

3.875 

10. 

95. 560 

4.  000 

135. 

516. 332 

3.925 

60. 

-314. 949 

3.  875 

10. 

105.823 

3.975 

110. 

526. 594 

3.  925 

60. 

-304. 686 

3.850 

-15. 

116.085 

3.  950 

85. 

536. 857 

3.900 

35. 

-294.  424 

3.  875 

10. 

126. 348 

3.950 

85. 

547. 120 

3.  925 

60. 

-284. 161 

3.875 

10. 

136.611 

3.950 

85. 

557. 383 

3.  900 

35. 

-273. 898 

3.  850 

-15. 

146.874 

3.925 

60. 

567. 645 

3.925 

60. 

-263. 635 

3.875 

10. 

157. 136 

3.950 

85. 

577. 908 

3.900 

35. 

-253.  373 

3.  875 

10. 

167.399 

3.925 

60. 

588. 171 

3.900 

35. 

-243. 110 

3.875 

10. 

177.662 

3.  925 

60. 

598. 433 

3.875 

10. 

-232. 847 

3.  850 

-15. 

187.924 

3.925 

60. 

608. 696 

3.875 

10. 

-222. 585 

3.  850 

-15. 

198. 187 

3.  925 

60. 

618. 959 

3.900 

35. 

-212. 322 

3.  850 

-15. 

208. 450 

3.925 

60. 

629. 222 

3.925 

60. 

-202. 059 

3.875 

10. 

218.713 

3.925 

60. 

639. 484 

3.875 

10. 

-191.796 

3.  875 

10. 

228. 975 

3.  925 

60. 

649. 747 

3.  900 

35. 

-181.534 

3.825 

-40. 

239. 238 

3.925 

60. 

660.010 

3.900 

35. 

-171.271 

3.  875 

10. 

249. 501 

3.925 

60. 

670. 272 

3.  900 

35. 

-161.008 

3.850 

-15. 

259. 763 

3.925 

60. 

680. 535 

3.850 

-15. 

-150.745 

3.  850 

-15. 

270. 026 

3.900 

35. 

690. 798 

3.  900 

35. 

-140.483 

3.850 

-15. 

280. 289 

3.875 

10. 

701.061 

3.875 

10. 

-130. 220 

3.  875 

10. 

290. 552 

3.900 

35. 

711.323 

3.  875 

10. 

-119.957 

3.875 

10. 

300.814 

3.925 

60. 

721.586 

3.875 

10. 

-109. 695 

3.  875 

10. 

311.077 

3.900 

35. 

731.849 

3.950 

85. 

-99. 432 

3.850 

-15. 

321.340 

3.875 

10. 

742. 112 

3.875 

10. 

-89. 169 

3.875 

10. 

331. 603 

3.  900 

35. 

752. 374 

3.875 

10. 

-78. 906 

3.875 

10. 

341.865 

3.875 

10. 

762. 637 

3.  850 

-15. 

-68. 644 

3.875 

10. 

352. 128 

3.900 

35. 

772. 900 

3.875 

10. 

-58.381 

3.850 

-15. 

362.391 

3.925 

60. 

783. 162 

3.875 

10. 

-48. 118 

3.  875 

10. 

372. 653 

3.  925 

60. 

793. 425 

3.900 

35. 

-37. 856 

3.850 

-15. 

382.916 

3.925 

60. 

803. 688 

3.900 

35. 

-27. 593 

3.850 

-15. 

393. 179 

3.900 

35. 

813.951 

3.875 

10. 

-17.330 

3.875 

10. 

403.442 

3.900 

35. 

824.213 

3.875 

10. 

*  -7.067 

3.925 

60. 

413.  704 

3.  875 

10. 

834. 476 

3.900 

35. 

*  3. 195 

4.025 

160. 

423.967 

3.925 

60. 

844. 739 

3.875 

10. 

*  13.458 

3.  900 

35. 

434. 230 

3.900 

35. 

855. 002 

3.  900 

35. 

23.721 

4.050 

185. 

444.493 

3.  925 

60. 

865. 264 

3.900 

35. 

33.  984 

4.  100 

235. 

454. 755 

3.900 

35. 

875. 527 

3.  875 

10. 

44. 246 

4.000 

135. 

465.018 

3.900 

35. 

885. 790 

3.850 

-15. 

54. 509 

4.050 

185. 

475.281 

3.  900 

35. 

Table  1.  Cervical  EP  Latency  Data  of  a  Selected  Peak:  Animal  „ 
AR  -  0761,  Experiment  0LX3OO9,  Acceleration  796  m/sec 
direction. 


•X 


f  Latencies  As  a  Function  of  Time,  Relative  to  Impact,  For  a  Selected 
(Experiment  #LX3009) . 


latency  value  In  Table  1  and  the  baseline  latency  value.  These 
shifts  (see  Table  1)  and  corresponding  times  were  then  analyzed 
using  the  FORTRAN  program  presented  in  the  Appendix. 


In  addition  to  the  data,  the  user  must  also  supply  the  program 
with  two  values  which  bracket  T  *,  where  T  is  the  exponential  time 
coefficient  in  model  (1).  A  lower  bound  of  -.1  sec  *  and  an  upper 
bound  of  0  sec  *  were  specified  for  T  *.  This  is  tantamount  to 
specifying  that  T  <  -10  seconds.  There  were  119  data  points,  of  which 
34  were  preimpact  and  85  were  postimpact. 

The  final  estimates  of  the  five  model  parameters  were  as  follows: 


A 


s  - 

-.04658660 

Usec/sec 

A 

B  - 

-8.25491618 

ysec 

A 

D  - 

66.28035173 

Usee 

/s 

A  - 

223.59247934 

Usee 

A 

T  - 

-65.27566084 

sec 

The  estimate  of  root  mean  square  deviation  was  19.017  Usee.  The 

2 

lack  of  fit  statistic,  MAR,  for  the  CCD  was  MAR  »  663.6  (Usee)  , 

A 

A 

which  was  1.6%  of  the  residual  sum  of  squares  [i.e.f  V(s»  b)  ] . 

A  plot  of  the  fitted  curve  is  shown  in  Figure  5. 

The  data  was  reanalyzed  using  the  seven  minute  time  interval 
starting  at  two  minutes  preimpact  and  ending  five  minutes  postimpact. 
There  were  39  data  points  in  this  interval,  of  which  11  were 
preimpact  and  28  were  postimpact.  The  estimated  parameters  were 


as  follows: 


5  ■  -.10731488  ysec/sec 

B  -  -6.45739421  usee 

6  -  78.29595060  usee 
A  -  205.91875792  usee 
T  -  -63.04820536  sec. 

The  estimate  of  root  mean  square  deviation  was  19.768  Usee. 

The  lack  of  fit  statistic,  MAR,  for  the  CCD  was  MAR  ■  372.6 
2 

(Usee)  ,  which  was  2.8%  of  the  residual  sum  of  squares. 
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APPENDIX 


11 

10 

6 

7 


8 


53 


IMPLICIT  REAL*8(A-H,  P-Z) 

REAL*8  XD  (600) ,  YD  (600) ,  T  (3) ,  XX  (600) ,  YY  (600) ,  Z  (6,  9) ,  T3LIM  (2) , 

&  XI  (9) ,  X2  (9) ,  YC  (9) ,  BETA  (6) ,  TPRE  (3) 

COMMON  DELT, El, E2, MAXIT,  MPRE,  MPOST,  XD,  YD,  T3LIM,  TOL 
DATA  Z/0.  0D0,  2.  0D0,  2.  0D0,  1.  0D0,  1.0D0,4.0D0, 

&  0.  0D0,  -2.  0D0,  2.  0D0,  1 . 0D0,  1 . 0D0,  -4.  0D0, 

&  0.  0D0,  -2.  0D0,  -2.  0D0,  1. 0D0,  1.  0D0,  4.  0D0, 

&  0.  0D0,  2.  0D0,  -2.  0D0,  1 . 0D0,  1 . 0D0,  -4.  0D0, 

&  16.  0D0,  0.  0D0,  0.  0D0,  -8.  0D0,  -8.  0D0,  0.  0D0, 

&  0. 0D0,  0.  0D0,  2. 0D0,  -1. 0D0, 3. 0D0, 0. 0D0, 

&  0.  0D0,  -2.  0D0,  0.  0D0,  3.  0D0,  -1.  0D0,  0.  0D0, 

&  0.  0D0,  2.  0D0,  0.  0D0,  3.  0D0,  -1 . 0D0,  0.  0D0, 

8,  0.  0D0,  0.  0D0,  -2.  0D0,  - 1 .  0D0,  3.  0D0,  0.  0D0/ 

Z (2, 7) =Z (2, 7) *DSQRT (2.  0D0) 

Z  (2, 8)  =Z (2, 8) *DSQRT (2.  0D0) 

Z (3, 6) =Z (3, 6) *DSQRT (2. 0D0) 

OPEN (6, FILE=’ PRN’ ) 

Z (3, 9)  =Z (3,  9) *DSQRT (2. 0D0) 

DO  10  1  =  1,6 
DO  11  J=l,  9 
Z  (I,  J)=Z(I,  J)/16.  0D0 
CONTINUE 
WRITE (*, 7) 

A  NONLINEAR  REGRESSION  PROCEDURE  FOR  EVOKED  POTENTIAL’, 
DATA  ANALYSIS  BY  CARL  A.  MAURO,  DESMATICS  INC.,’, 

PO  BOX  618,  STATE  COLLEGE,  PA  16804  (814-238-9621).’, 
BASED  ON  DESMATICS  TECHNICAL  REPORT  NO.  112-18’, 

PROGRAM  FITS  A  REGRESSION  FUNCTION  OF  THE  FORM’, 
Y=B+St+h(t) CD+Aexp (t/T) 3  WHERE:’ , 
t  =  TIME  RELATIVE  TO  IMPACT (SECS) ’ , 

Y  =  LATENCY  SHIFT (MICROSECS)  RELATIVE  TO  BASELINE  AEP* , 
h (t )  =  ZERO  FOR  t  <0  AND  UNITY  FOR  t>=0*, 


FORMAT ( 

/ 
/ 

// 
/ 
/ 

// 
/ 
/ 
/ 
/ 
/ 
/ 
/ 

PAUSE 
WRITE (* 
FORMAT ( 


/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

//// 


B  =  AMPLITUTE  INSTATANEOUSLY  PRIOR  TO  IMPACT (MICROSECS) ’ 
S  =  SLOPE  OF  LINEAR  BASELINE  SHIFT (MICROSECS/SEC) ’ , 

D  =  SHIFT  INDUCED  BY  IMPACT (MICROSECS) ’ , 

A  =  EXPONENTIAL  AMPLITUDE (MICROSECS) ’ , 

T  =  EXPONENTIAL  TIME  COEFFICIENT (SECS) ’////// ) 


8) 

///// 

NOTE: 


(1) 


(2) 


(3) 


(4) 


DATA 

WITH 

EACH 


IS  READ  FREE-FORMAT  FROM  THE  FILE*, 
FILENAME  DATA. EP’ , 

RECORD  IN  DATA. EP  MUST  CONTAIN  THE  TIME*, 
IN  SECONDS  OF  THE  DATA  POINT,  FOLLOWED  BY  THE’ 
CORRESPONDING  Y  VALUE  IN  MICROSECONDS’, 

DATA  ARRAYS  ARE  DIMENSIONED  TO  ACCOMODATE’, 

A  MAXIMIMUM  OF  600  DATA  POINTS’, 

IN  ADDITION  TO  THE  DATA,  USER  MUST  SUPPLY  TWO’ 
VALUES  WHICH  BRACKET  THE  RECIPROCAL  OF  T»  , 
PLEASE  INPUT  LOWER  LIMIT  FOR  1/T  (IN  1/SECS)’) 

READ(*,  *)  T3LIM ( 1 ) 

WRITE (*,  53) 

FORMAT (//’  PLEASE  INPUT  UPPER  LIMIT  FOR  1/T  (IN  1/SECS)’) 

READ (*, *)  T3LIM (2) 

IF  (T3LIM (2) . LE. T3LIM ( 1 ) )  GO  TO  6 
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*  -  *  .  *  . 


SF*10. 0D0 
TOL*.  0001 
DELT=2.**(-8) 

El*. 0001 
E2*. 001 
MAX IT=25 
MPRE=0 
MP0ST=0 

C  >>>>>  INPUT  DATA 

WRITE  <*, 54) 

54  FORMAT (/’  READING  DATA  .....’) 

OPEN ( 1 , FILE*’  DATA.  EP’  ) 

DO  5  I*  1,601 

READ <1,*, END* 15)  XD(I),YD(1) 

5  CONTINUE 

15  M=I-1 

CLOSE  (1) 

C  >>>>>  SORT  DATA 

ND*M 
NM=M— 1 

DO  20  1*1, NM 
ND*ND-1 
DO  21  J*1 , ND 

IF  (XD(J).LE. XD(J+1> )  GOTO  21 
A*XD(J) 

B*YD(J) 

XD ( J) *XD  (J+l ) 

YD ( J) *YD  ( J+l ) 

XD  < J+l ) =A 
YD( J+l ) =B 
21  CONTINUE 
20  CONTINUE 

C  >>>>>  SEPARATE  PRE  FROM  POST  IMPACT  DATA 
DO  30  1=1, M 

IF  (XD(I).  BE.0.  )  GO  TO  31 

MPRE=MPRE+1 

XX (MPRE) =XD ( I ) 

YY (MPRE) =YD ( I ) 

GO  TO  30 

31  MPOST =MPOST  + 1 

30  CONTINUE 

C  >>>>>  CHECK  FOR  SUFFICIENT  DATA  POINTS 

IF  (MPRE.LT. 3.0R. MPOST.LT. 6)  STOP  ’INSUFFICIENT  DATA’ 

WRITE (*,55)  MPRE, MPOST 

55  FORMAT (/’  NUMBER  OF  PREIMPACT  DATA  POINTS*’ , 14, 

ft  /’  NUMBER  OF  POSTIMPACT  DATA  POINTS*’ , 14, /’  COMPUTING  . ’) 

C  >>>>>  COMPUTE  INITIAL  ESTIMATES  OF  BASELINE  SLOPE  AND  AMPLITUDE ( I NT) 
CALL  SREG (MPRE,  XX,  YY, BZ, B1 ) 

C  >>>>>  COMPUTE  ESTIMATED  STANDARD  ERRORS  OF  INITIAL  BZ  AND 
Bl (SDBZ, SDB1 ) 

PSSE*0. 0D0 
XXBAR*0. 0D0 
DO  35  1*1, MPRE 

XXBAR*XXBAR+XX ( I ) 

35  PSSE*PSSE+ ( YY ( I ) -BZ-B1*XX ( I ) ) *#2 


i  kSaWfckjio  iw 


/.  -v 


-V- W 
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n  n 


XXBAR=XXBAR/MPRE 
SD=DSQRT ( PSSE/ (MPRE-2.  DO) ) 

SXX=0. ODO 
DO  36  1=1, MPRE 

36  SXX=SXX+ (XX ( I ) -XXBAR) ##2 

A=1 . / ( 1 . DO*MPRE )  +  ( X  XBAR**2) /SX  X 
SDB1=SD/DSQRT (SXX) 

SDBZ=SD*DSQRT (A) 

>>>>>  START  CENTRAL  COMPOSITE  DESIBN 
>>>>>  DEFINE  DESIGN  MATRIX 
40  XI (1)=B1+SDB1/SF 
XI (2)=B1-SDB1/SF 
XI (3) =X1 (2) 

XI (4) =X1  (1) 

X 1 (5) =B1 
XI (6) =B1 

X 1 ( 7 ) =B1 -DSQRT ( 2. ODO ) *SDB 1 / SF 
XI (8) =B1+DSQRT (2. ODO) *SDB1/SF 
XI (9)=B1 

X2(1)=BZ+SDBZ/SF 
X2(2) =X2 ( 1 ) 

X2 ( 3  > “BZ-SDBZ /SF 
X2  (4)  =X2  (3) 

X2 (5) =BZ 

X2(6) =BZ+DSQRT (2. ODO) #SDBZ/SF 

X2(7)=BZ 

X2(8) =BZ 

X2 (9) =BZ-DSQRT (2. ODO) *SDBZ/SF 
DO  100  1=1,9 

CALL  EGAD (XI (I),  X2 ( I ) ,  T,  YC ( I ) ) 

IF  (I.NE.5)  GOTO  100 
TPRE  < 1 ) =T ( 1 ) 

TPRE  <2) =T (2) 

TPRE (3) =1. ODO/T (3) 

100  CONTINUE 

C  >>>>>  DETERMINE  BETAS  FOR  QUADRATIC  FIT 
DO  105  1  =  1,6 
BETA ( I ) =0.  ODO 
DO  110  J=1 , 9 

HO  BETA  ( I )  =BETA  ( I )  +Z  <  I,  J)  *YC  ( J) 

105  CONTINUE 

DETER=BET A ( 4 >  *BETA ( 5 ) -BETA (6) *BETA(6) /4. ODO 
C  >>>)>  DETERMINE  STATIONARY  POINT 

STAT 1 -BETA (2) *BETA (5) -BETA (3) *BETA (6) /2. ODO 
STAT1—STAT1/  (2.  ODO*DETER) 

STATZ-BETA (3) *BETA (4) -BETA  <2) #BETA (6) /2. ODO 
STATZ— STATZ/  (2.  ODO*DETER) 

MDB=DSQRT< (STATZ**2+STAT1#*2) /2. ODO) 

C  >>>>>  TRANSFORM  TO  ORIGINAL  SPACE 

STAT 1 «B 1 +STAT 1 *SDB 1 /SF 
STATZ=BZ+STATZ#SDBZ/SF 
CALL  EGAD ( STAT 1, STATZ, T,FSTAT) 

RMS-DSQRT (YC (5) / <M-5. ODO) ) 

WRITE  (6,150)  Bl,  BZ, TPRE ( 1 ) , TPRE (2) , TPRE (3) , YC (5) , RMS 
ISO  FORMAT (/////’ 1CENTER  POINT  OF  CENTRAL  COMPOSITE  DESIGN 


&  /*  S=*,F15.8,/’  B=* ,  F15.  8,  /’  CONDITIONAL  ESTIMATES?* , 

&  /»  D=*,F15.8,  /»  A»’,F15.  8, /’  T=’,F15.8, 
ft  /*  SUM  OF  SQUARED  RESIDUALS'3* ,  FI 5.  8, 
ft  /»  ROOT  MEAN  SQUARE-’ , F15. 7) 

RMS-DSQRT (FSTAT/ <M-5.  0D0) ) 

X MAR-0. 0D0 
DO  153  1  =  1,9 

YCP=BETA<1)+BETA<2)#X1  < I ) +BETA (3) *X2 ( I ) +BETA <4) *X1  <I)#X1 <I) 
ft  +BETA (5) *X2 ( I ) *X2 ( I ) +BETA  <6) *X1 < I ) #X2 ( I ) 

XMAR-XMAR+DABS  < YC  < I ) -YCP) 

153  CONTINUE 

XMAR-XMAR/9.  0D0 

WRITE  <6,160)  STAT1,STATZ,T<1),T<2),  1.0D0/T<3>, FSTAT,  RMS,  XMAR 
160  FORMAT (/////*  FINAL  ESTIMATES  OF  MODEL  PARAMETERS?’, 
ft  /’  S=’,F15.8, /’  B=’,F15.8, /’  D=* ,  F15. 8,  /’  A=’,F15.8, 
ft  /’  T=* , F15. 8, /*  SUM  OF  SQUARED  RESIDUALS-* , F15. 6, 
ft  /»  ROOT  MEAN  SQUARE-’ , F 15. 7, 

ft  /’  LACK  OF  FIT  STATISTIC  FOR  CENTRAL  COMPOSITE  DESIGN i* , 
ft  ’  MAR-*  ,F15.6//> 

IF  (UDB. LE. DSQRT  <2. 0D0) )  GO  TO  200 
WRITE <6, 59) 

59  FORMAT <’  NOTE:  STATIONARY  POINT  IS  OUTSIDE  BOUNDARIES  OF*, 

ft  /*  CENTRAL  COMPOSITE  DESIGN.  DESIGN  WILL  BE  RECENTERED*, 

ft  /’  AND  NEW  RESPONSE  SURFACE  FIT  TO  (S,B)  SPACE.*, 

ft  /’  COMPUTING . ’  ) 

B1-STAT1 
BZ-STATZ 
GO  TO  40 
200  STOP 
END 

SUBROUTINE  EGAD (PSLP,  PINT,  T,  FUNV) 

IMPLICIT  REAL*8<A-H, P-Z) 

REAL*8  X  <600) ,  Y <600) ,  T  <3) ,  TH<3) , YR  <600) ,  JM<600,  3) , TT <3) , 
ft  JR<3),G<3,3),D<3),  Z<3,3),  WK<3),  V  <3,  3,  3) ,  GMINV  <3,  3), 
ft  XSL  <600) , XD<600) , YD<600) ,  T3LIM  <2) 

COMMON  DELT, El, E2, MAX  IT, MPRE, MPOST, XD, YD, T3LIM, TOL 
C  >>))>  FORM  ADJUSTED  POST  IMPACT  DATA  ARRAYS 

DO  4  1=1, MPOST 
X  <1 )=XD< I+MPRE) 

Y  < I ) -YD  < I+MPRE) -PINT-PSLP*XD  < I+MPRE) 

4  CONTINUE 

C  >))))  COMPUTE  STARTING  ESTIMATES  OF  T < 1 > , T (2) , T <3) 

A-T3LIM  < 1 ) 

B— T3LIM<2) 

XMIN  *  A 
IER  -  129 

IF  <B  .LE.  A)  GO  TO  24 
IER  -  130 

IF  <TOL  .GE.  <B-A) )  GO  TO  24 
IER  «  0 

C  =  <3. 0D0-DSQRT <5. 0D0) ) /2. 0D0 
H  -  C*  <B-A) 

VI  -  A+H 
V2  -  B-H 
FA  -  F  <A,  X,  Y> 
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FB  =  F  (B,  X,  Y) 

FV1  -  F(V1,X,Y> 

FV2  -  F(V2,  X,Y) 

S  CONTINUE 

IF  (A  .BE.  VI  .OR.  VI  .BE.  V2  .OR.  V2  .BE.  B>  BO  TO  23 
IF  (FV1  .BE.  FV2)  60  TO  10 
IF  (FV2  .ST.  FB)  BO  TO  25 
B  -  V2 

IF  (TOL  .BE.  (B-A) )  SO  TO  15 
FB  «  FV2 
V2  -  VI 
FV2  «  FV1 
H  «  C* (B— A) 

VI  *  A+H 

FV1  -  F (VI, X, Y) 

SO  TO  5 

10  IF  (FV1  .ST.  FA)  SO  TO  21 
A  -  VI 

IF  (TOL  .BE.  (B-A))  SO  TO  20 
FA  -  FV1 
VI  -  V2 
FV1  -  FV2 
H  -  C*(B-A) 

V2  »  B-H 

FV2  «  F(V2,  X,Y) 

SO  TO  5 
15  XMIN  -  VI 

IF  (FA  .LT.  FV1 )  XMIN  -  A 
BO  TO  27 

20  XMIN  -  V2 

IF  (FB  .LT.  FV2)  XMIN  *  B 
BO  TO  27 
25  XMIN  »  V2 
A  *  VI 
BO  TO  22 

21  XMIN  -  VI 
B  3  V2 

22  IER  3  131 
BO  TO  24 

23  IER  ■  132 
XMIN  -  A 

IF  (FB  .LT.  FA)  XMIN  «  B 
24  CONTINUE 

IF  (IER  .EQ.  0)  SO  TO  27 
WRITE  (*,26)  IER 

26  FORMATS  TERMINAL  ERROR  FROM  BOLDSECTION.  ERROR  CODE-  ’,13) 
STOP 

27  CONTINUE 
T (3) -XMIN 
DO  31  1-1, MPOST 
31  XSL ( I ) »DEXP(T (3) *X ( I )  ) 

CALL  SREB(MPOST, XSL,  Y,  BINT,  BSLP) 

T ( 1 ) -BINT 
T(2)-BSLP 

C  >>>>>  START  ITERATION  LOOP 


•“  ,*  •  ^  ,•  »  -  «  *  .»  •  ««•»*<  •  -  -  *i*i  *•  •  » 
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DO  100  I IT-1, MAXIT 

>>>>>  COMPUTE  RESIDUAL  VECTOR  (YR)  AND  JACOBIAN  MATRIX  <JM> 
SS-0. 

DO  30  J— 1,  MPOST 

YR ( J) =Y ( J) -T ( 1 ) -T (E) *DEXP (T  <3>  #X  ( J)  ) 

SS-SS+YR(J)*YR<J> 

JM  ( J,  1>=-1. 

JM(  J,  2) =-DEXP  (T  (3)  *X  <  J) > 

JM(J, 3)— X(J)*T(2)*DEXP<T<3)*X(J) > 

30  CONTINUE 

>>>>>  FORM  MATRIX  PRODUCT  JR»(TRAN  OF  JM)*YR 
JR  ( 1 )  =0. 

JR  (2)  -0. 

JR (3) =0. 

DO  35  J-l, MPOST 
JR ( 1 ) =JR  < 1 ) -YR ( J) 

JR (2) -JR (2) +JM ( J,  2) *YR ( J) 

JR (3) -JR (3) +JM ( J, 3)#YR<J> 

35  CONTINUE 

>>>>>  COMPUTE  HESSIAN  MATRIX  (G) 

DO  40  1-1,3 
DO  45  J-l,  3 
0<I,  J>=0. 

DO  50  L-l, MPOST 
G<I,J)=G(I,J)  +JM  (L,  I)*JM(L, J) 

IF  (I+J.EQ.  6)G(I, J)-G<I,  J)-YR(L)#X<L)*X(L>*T(2)*DEXP(T(3)*X(L) ) 
IF  (I+J.EQ.5)  G(I,J)=G(I,J)  -YR  <L)  *X  <L>  *DEXP (T  <3) *X (L)  ) 

50  CONTINUE 
45  CONTINUE 
40  CONTINUE 

>>>>>  CALCULATE  INVERSE  OF  MODIFIED  HESSIAN  MATRIX 
CALL  EIGRS  (G,  3,  1 1,  D,  Z,  3,  UK,  IER) 

DO  60  J=l,  3 

60  D ( J) -DMAX1 ( ABS (D  ( J)  ) , DELT) 

DO  62  1-1,3 
DO  63  J-l,  3 
DO  64  K-l,  3 

64  V  <  I ,  J,  K)  -Z  <  J,  I )  *Z  (K,  I  > 

63  CONTINUE 

62  CONTINUE 

DO  65  J-l, 3 
DO  66  K= 1 , 3 
GMINV  ( J,  K)  =0. 

DO  67  1-1,3 

67  GMINV <J,K> -GMINV (J,K)+V( I,  J,  K) /D(I> 

66  CONTINUE 

65  CONTINUE 

>>>>>  UPDATE  ESTIMATES 
DO  70  1-1,3 
TH  < I >  -0. 

DO  75  J-l,  3 

75  TH(I)»TH<I>+GMINV<I, J)*JR(J) 

70  CONTINUE 

>>>>>  STEP  SIZE  CALCULATION 

ALP-2. 
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1=1 

l  ALP=ALP/2. 

DO  81  J=l, 3 
TT(J>=T<J)-ALP*TH(J> 

CONTINUE 

SST=0. 

DO  82  J=l, MPOST 

S1=Y<J)-TT<1)-TT<2)*DEXP<TT(3>*X<J> ) 

SST =SST +S 1 *S 1 
>  CONTINUE 

IF  (SST.LT.SS)  GOTO  85 
1=1+1 

IF (I.  EQ. 10)  GOTO  85 
GOTO  80 
5  CONTINUE 

>>>>>  CHECK  FOR  TERMINATION 
ITERM-1 
DO  90  1  =  1,3 
F 1 = ABS ( T ( I ) -TT  < I ) ) 

F2=E 1  * ( ABS  <  T ( I ) ) +E2 ) 

IF  (Fl.GT. F2)  ITERM=0 
T  ( I )=TT ( I ) 

9  CONTINUE 

IF  (I TERM. EQ. 1)  GOTO  101 

90  CONTINUE 

91  IF ( ITERM. EQ.  1 )  NITT-IIT 

IF ( ITERM.  EQ.  0)  STOP  ’MAX  ITERATIONS  EXCEEDED  IN  POST  MODULE’ 
FUNV=0. 

DO  150  1=1 , MPRE 

FUN V=FUN V+  <  YD  < I ) - P I NT- PSL  P*XD< I) >**2 
50  CONTINUE 

DO  160  J=l, MPOST 
I=J+MPRE 

A=PINT+PSLP*XD(I) 

B=T ( 1 ) +T (2) *DEXP (T (3) *XD  < I ) ) 

FUNV=FUNV+ ( YD ( I ) -A-B ) **2 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  SREG (N, X, Y, BINT, BSLP) 

IMPLICIT  REAL*8(A-H, P-Z) 

REAL*8  X (600) ,  Y (600) 

XBAR=0. 

YBAR=0. 

DO  5  1  =  1,  N 
XBAR=XBAR+X ( I ) 

YBAR=YBAR+Y ( I ) 

CONTINUE 

XBAR-XBAR/N 

YBAR=YBAR/N 

A=0. 

B»0. 

DO  10  1-1,  N 
C=X (I)-XBAR 
A=A+Y ( I ) *C 
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B=B+C**2 
10  CONTINUE 
BSLP=A/B 

BINT=YBAR-BSLP*XBAR 

RETURN 

END 

SUBROUTINE  EIGRS  (A,  N, JOBN, D,  Z, IZ, WK, IER) 

DOUBLE  PRECISION  A (9) ,  D (3) , WK (3) , Z <3,  3) 

INTEGER  N,  JOBN,  IZ,  IER,  IJ,  JI 

INTEGER  JER,  NA,ND,  IIZ,  I,  J,K 

DOUBLE  PRECISION  TEN,  ZERO,  ONE,  THOUS 

DATA  ZERO,  ONE/0. 0D0,  1. 0D0/, TEN/ 10. 0D0/,  THOUS/1000.  0D0/ 

IER  =  0 
JER  =  0 
K  =  1 
JI  »  N-l 
IJ  -  1 
DO  10  J=1,N 
DO  5  1=1, J 

A<K)  =  A ( I J) 

IJ  ■  IJ+1 
K  =  K+l 
5  CONTINUE 

IJ  =  IJ  +  JI 
JI  *  JI  -  1 
10  CONTINUE 

NA  =  (N*(N+1) )/2 
ND  =  1 

CALL  EHOUSS  (A, N, D, WK, WK) 

IIZ  =  IZ 
DO  55  1=1, N 
DO  50  J=1,N 

Z(I,J>  -  ZERO 
50  CONTINUE 

Z(I, I)  -  ONE 
55  CONTINUE 

CALL  EQRT2S  (D, WK, N, Z, I I Z, JER) 

IF  (JER. GT. 128)  GO  TO  9000 
CALL  EHOBKS  (A,  N,  1,  N,  Z,  I Z ) 

9000  CONTINUE 

IF  (JER.  EQ.0)  GO  TO  9005 

STOP  *  NOT  ALL  EIGENVALUES  COULD  BE  COMPUTED* 

9005  RETURN 
END 

SUBROUTINE  EHOBKS  (A, N,  Ml , M2, Z, IZ ) 

DIMENSION  A (9) , Z (3, 3) 

DOUBLE  PRECISION  A, Z,H, S 
DO  25  1=2,  N 
L  -  1-1 
IA  -  (I*L)/2 
H  =  A ( I A+I ) 

IF  (H.EQ.0.  D0)  GO  TO  25 
DO  20  J  =  Ml, M2 
S  =  0. 0D0 
DO  10  K  >  1, L 
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S  =  S+A ( IA+K) *Z  <K, J> 

10  CONTINUE 

S  -  S/H 
DO  15  K=1,L 

Z(K,J)  »  Z (K,  J> -S*A ( IA+K) 
15  CONTINUE 

20  CONTINUE 
25  CONTINUE 
RETURN 
END 

SUBROUTINE  EHOUSS  (A, N, D, E,  E2) 


DIMENSION  A (9) ,  D (3) , E (3 

DOUBLE  PRECISION  A, D, E, E2f ZERO 
DATA  ZERO/0. 0D0/, 0 

NP1  -  N+l 
NN  *  (N*NPl)/2-l 
NBEG  =  NN+l-N 
DO  70  II  *  1, N 
I  =  NP1-II 
L  =  1-1 
H  =  ZERO 
SCALE  ■  ZERO 
IF  (L  .LT.  1)  GO  TO  10 
NK  -  NN 
DO  5  K  -  1,L 

SCALE  =  SCALE+DABS (A (NK) > 
NK  *  NK— 1 
CONTINUE 

IF  (SCALE  .NE.  ZERO)  GO  TO  15 
E(I)  -  ZERO 
E2(I)  »  ZERO 
GO  TO  65 
NK  =  NN 

SCALE 1  «  ONE /SCALE 
DO  20  K  *  1,L 

A (NK)  *  A (NK) *SCALE1 
H  =  H+A (NK) *A (NK) 

NK  *  NK-1 
CONTINUE 

E2(I)  *  SCALE*SCALE*H 
F  =  A  (NN) 

G  =  -DSIGN (DSQRT (H) , F) 

E(I)  =  SCALE*G 
H  ®  H-F*G 
A  (NN)  =»  F-G 
IF  (L  .EQ.  1)  GO  TO  55 
F  ■  ZERO 
JK1  -  1 
DO  40  J  «  1, L 
G  -  ZERO 
IK  -  NBEG+1 
JK  »  JK1 
DO  25  K  =  1,J 

G  -  G+A(JK)#A(IK) 


A(9),D(3),E(3),E2(3) 

A,  D,  E,  E2f  ZERO,  H,  SCALE, ONE, SCALE1, F 
ZERO/0. 0D0/, ONE/1. 0D0/ 


IK  *  IK+1 
25  CONTINUE 

JP1  =  J+l 

IF  (L  .LT.  JP1)  SO  TO  35 
JK  -  JK+J-1 
DO  30  K  -  JP1,  L 

6  *  G+A(JK)*A(IK) 

JK  =  JK+K 
IK  =  IK+1 
30  CONTINUE 

35  E  (J)  *  S/H 

F  «  F+E ( J>  *A  <NBEG+J) 

JK1  =  JK1+J 
40  CONTINUE 

HH  =  F/<H+H) 

JK  =  1 

DO  50  J  *  1,L 

F  =  A (NBES+J) 

G  =  E(J)-HH*F 
E  ( J)  =  G 
DO  45  K  -  1,J 

R<JK>  »  R  ( JK) -F*E (K) -G#A (NBEG+K) 

JK  =  JK+1 
45  CONTINUE 

50  CONTINUE 

55  DO  60  K  ■  1,L 

fl < NBEG+K)  *  SCRLE*fl (NBEG+K) 

60  CONTINUE 

65  D(I)  -  fl (NBEG+I ) 

A(NBEG+I )  -  H*SCRLE*SCRLE 
NBEG  ■  NBEG-I+1 
NN  -  NN-I 
70  CONTINUE 
RETURN 
END 

SUBROUTINE  EQRT2S  <D, E, N, Z, I Z, IER) 

DIMENSION  D (3) , E (3) , Z (3, 3) 

DOUBLE  PRECISION  D,  E,  Z,  B,  C,  F,  G,  H,  Py  R,  S,  RDELP,  ONE,  ZERO 
DATA  ZERO, ONE/0. 0D0, 1 . 0D0/ 

IER  -  0 
RDELP»a.**<-5a> 

DO  5  1=2, N 
E (1-1 ) =E ( I ) 

5  CONTINUE 

E (N)  =  ZERO 
B  -  ZERO 
F  »  ZERO 
DO  60  L-1,N 

J  =  0 

H  =  RDELP* (DABS (D (L) ) +DABS (E (L) ) ) 

IF  (B.  LT.  H)  B  =  H 
DO  10  M=L,N 
K=M 

IF  (DABS(£(K) )  .LE.  B)  GO  TO  15 
10  CONTINUE 


15 

M  a 

■  K 

IF 

(M.  EQ.L) 

GO 

TO 

55 

20 

IF 

(J  .EQ. 

30) 

GO 

TO 

J  *  J+l 

LI 

*  L+l 

S  =  D (L) 

P  *  (D(L1)-G>/<E<L)+E<L) > 

R  *  DSQRT (P*P+ONE) 

D(L)  -  E  <L)  / (P+DSIGN  <R,  P) > 

H  »  G-D(L) 

DO  25  I  »  LI, N 
D ( I )  -  D  < I ) -H 
25  CONTINUE 
F  *  F+H 
P  *  D(M) 

C  »  ONE 
S  *  ZERO 
MM1  *  M-l 
MMI  PL  =  MM1+L 
IF  (L.GT.MM1)  GO  TO  50 
DO  45  II=L, MMI 
I  *  MMI PL— 1 1 
G  =  C*E ( I ) 

H  a  C#P 

IF  (DABS(P).LT.  DABS(E(I>)>  GO  TO  30 
C  *  E(I)/P 
R  *  DSQRT <C*C+ONE> 

E  ( 1+1 )  *  S*P*R 
S  *  C/R 
C  =  ONE/R 
GO  TO  35 

30  C  «  P/E (I) 

R  =  DSQRT <C*C+QNE> 

E  ( 1  +  1 )  »  S*E ( I >  *R 
S  *  ONE/R 
C  «  C*S 

35  P  «  C*D(I)-S*G 

D(I+1)  *  H+S*(C*G+S#D(I) ) 

IF  (IZ  .LT.  N)  GO  TO  45 
DO  40  K-1,N 

H  *  Z  (K, I+l) 

Z(K, 1+1)  »  S*Z(K, I)+C#H 
Z(K,  I)  ■  C*Z (K, I > -S*H 
40  CONTINUE 

45  CONTINUE 

50  E(L)  »  S*P 

D(L)  *  C*P 

IF  (DABS  (E  (L) )  .GT.B)  GO  TO  20 
55  D(L)  *  D (L)  +  F 

60  CONTINUE 

DO  80  1*1,  N 

K  *  I 
P  ■  D  ( I ) 

IP1  -  1+1 

IF  (IPl.GT.N)  GO  TO  70 
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DO  65  J»IP1,N 

IF  <D(J)  .BE.  P)  BO  TO  65 
K  -  J 
P  »  D  < J) 

65  CONTINUE 

70  IF  (K.EQ.  I)  BO  TO  80 

D<K)  -  D ( I ) 

D  (I )  -  P 

IF  (IZ  .LT.  N)  SO  TO  80 
DO  75  J  »  1,N 
P  «  Z(JfI) 

Z<J,I>  -  Z  ( J,  K) 

Z(J,K)  -  P 
75  CONTINUE 
80  CONTINUE 
BO  TO  9005 
85  IER  =  128+L 
9000  CONTINUE 
9005  RETURN 
END 

DOUBLE  PRECISION  FUNCTION  F (A, X, Y) 

IMPLICIT  REAL*8<A-H, P-Z) 

REAL*8  X(600), Y(600), XX (600) 

COMMON  DELTf El , E2, MAX  IT, MPRE, MPOST, XD, YD, T3LIM, TOL 
YBAR=0. 0D0 
DO  5  1=1, MPOST 
YBAR«YBAR+Y <1 ) 

5  XX  < I ) =DEXP (A*X ( I) ) 

F-0. 0D0 

IF  (A.EO.0.0D0)  BO  TO  100 
CALL  SRE6< MPOST, XX, Y, BINT,  BSLP) 

DO  10  I«l, MPOST 
F=F+< Y  < I )— BINT-BSLP*XX (I) )**2 
10  CONTINUE 
RETURN 

100  YBAR*YBAR/MPOST 
DO  20  I -1, MPOST 
20  F=F+ ( Y ( I ) - YBAR ) **2 
RETURN 
END 
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